Global phase diagram of three-dimensional extended Boson Hubbard model - a 
continuous time Quantum Monte Carlo study 



O 

(N 

< 

in 

(N 



C/2 



o 
o 



> 

00 

o 



X 



Bin Xi^'^, Fei Ye^, Weiqiang Chen^, Fuchun Zhang^ and Gang Su^ 
of Physical Sciences, Graduate University of Chinese Academy of Sciences, Beijing 100049, China 
^ College of Material Sciences and Optoelectric Technology, 
Graduate University of Chinese Academy of Sciences, Beijing 100049, China 
Centre of Theoretical and Computational Physics and Department of Physics, 
the University of Hong Kong, Hong Kong, China and 
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China 

(Dated; January 20, 2013) 

We present the global phase diagram of the extended boson Hubbard model on a simple cubic 
lattice by quantum Monte Carlo simulation with worm update algorithm. Four kinds of phases are 
supported by this model, including superfluid, supersolid, Mott, and charge density wave (CDW) 
states, which are identified in the phase diagram of chemical potential /i versus nearest neighbor 
interaction V . By changing the chemical potential, a continuous transition is found from the Mott 
phase to a superfluid phase without breaking the translational symmetry. For an insulating CDW 
state, adding particles to it gives rise to a continuous transition to a supersolid phase, while removing 
particles usually leads to a first-order one to either supersolid or superfluid phase. By tuning the 
nearest neighbor interaction, one can realize the transition between two insulating phases, Mott and 
CDW with the same particle density, which turns out to be of the first-order. We also demonstrate 
that a supersolid phase with average particle density less than 1/2 can exist in a small region of 

— V phase diagram. 

PACS numbers: 67.80.kb, 75.40.Mg, 02.60.-x, 64.60.-i 



I. INTRODUCTION 

Lattice models of interacting bosons and fermions such 
as Hubbard model and its various generalizations are usu- 
ally strongly correlated systems exhibiting various phases 
with competing orders, which are of fundamental inter- 
est in fields of both condensed matter and cold atomic 
physics. Interests on both types of Hubbard models 
are renewed recently, since they can be realized in cold 
atomic gases loaded in optical lattices (for a review see 
Refs. [l| and and references therein). Unlike fermions, 
there is a natural superfluid order for free bosons at zero 
temperature driven by the kinetic energy. When the in- 
teraction is switched on, the bosons are likely to be lo- 
calized in various crystalline patterns, which may coexist 
with superfluid order to give a realization of intrigu- 
ing "supersolid" state that has been pursued for decades 
since 1950s [7l-[To|. Recently, people have observed the 
non-classical rotational inertia in solidified ''He [ill, [13 
implying a possible supersolid state, which, in spite of the 
controversy over this topic, also triggers extensive studies 
on various boson Hubbard models. 

Experimentally, the boson Hubbard model can be 
used to mimic the granular superconductors, where the 
Cooper pairs are described as bosons, which has been 
studied by Fisher et al. [ll] two decades ago, where 
with only on-site repulsive interaction they showed that 
bosons can form either Mott insulating state with integer 



filling or superfluid state. Recent experimental progress 
in cold atomic system provides another realization of bo- 
son Hubbard model by loading atoms into an optical lat- 
tice with possible long range interactions through dipole 
interaction (l4l - [l6ll, o r mediated by other intermediate 
states or fermions [17l - ll9| . In addition, the boson models 
also share similarities with quantum magnets, e.g., the 
uniaxial magnetization corresponds to insulating states 
of boson Hubbard model (e.g. Ref. [20l ). while the easy- 
plane magnetization corresponds to the superfluid state. 
Hence, the studies on the boson Hubbard model may 
shed light on some common issues of strongly correlated 
lattice models. 

Generally speaking, boson models with interactions 
at zero temperature have two principal phases: (i) the 
superfluid and (ii) the incompressible insulating state, 
which are favored respectively by kinetic and interaction 
energies, and can coexist. Depending on the features of 
interaction terms, there are several types of insulating 
phases, such as Mott, valence bond crystal, and charge 
density wave (CDW). Note that we in this article define 
the incompressible states with oscillating density profile 
as CDW, though the bosons may not carry charges. 

The extended boson Hubbard (EBH) model with on- 
site (U) and nearest neighbor (V) interactions is a mini- 
mal model in favor of CDW and supersolid phases, which 
has the form of 
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where {hi) is the creation (annihilation) bosonic op- 
erator at site i, t is the hopping amphtude, fii = blbi 
is the particle number, fi is the chemical potential, and 
runs over all nearest neighbors. Recently, Hamil- 
tonian Eq. ^ and its hard-core version (equivalent to 
the quantum spin-1/2 XXZ model) with different under- 
lying lattices have been extensively studied in different 
parameter regimes [sl IgI. [2ll - l28i] . However, a global phase 
diagram of the three-dimensional (3D) EBH model [Eq. 
([T])] is still of lack. As there is no sign problem for the 
EBH model, the quantum Monte Carlo (QMC) simula- 
tion is the most convenient tool for this purpose. The 
worm algorithm (29l - l3l| will be invoked to study Hamil- 
tonian ([T]) on a simple cubic lattice, together with other 
perturbation and mean-field approaches. 

The system described by the EBH model can give rise 
to a charge ordered crystal at commensurate fillings. The 
first one is for half filling p = 0.5, and the correspond- 
ing solid state is labeled as CDW I. Doping particles into 
this state can lead to a supersolid state [a, However, 
as shown in Ref. doping holes into it acts quite dif- 
ferently, which may not result in a supersolid state with 
p < 0.5, but a phase separation between superfluid and 
CDW I states, which signals a first-order phase transi- 
tion. Their argument is based upon the following two ob- 
servations. (I) Taking one particle out of a perfect CDW 
crystal with half filling costs almost no potential energy, 
but only chemical potential. At the same time, the hop- 
ping hole also gains a kinetic energy which is quadratic 
in t (~ t^). For a perfect CDW crystal, these three pro- 
cesses are balanced, so one cannot take one particle out. 
(H) The CDW phase breaks the translational symme- 
try, leading to a two-fold degenerate ground state. If 
holes are doped into the domain wall between these two 
degenerate phases, the kinetic energy gained is propor- 
tional to t. Hence, the CDW phase is unstable toward 
the domain wall formation if the hole density exceeds 
for L'^ lattice, though it is still stable against los- 
ing one particle. This argument perfectly explains the 
first-order phase transition from the CDW I to super- 
fluid state with p < 0.5, but it fails in two circumstances. 
The first is that in one dimension the kinetic energy is 
always linear in and the corresponding transition is of 
the Kosterlitz-Thouless type j6|] . The other is that if V is 
comparable to t the kinetic energy of holes is also linear 
in t, which may result in the supersolid phase with the 
particle density less than half filling (see Sec. IIIBt . This 
can be verified by the mean- field calculations |26l fs^ . 

At unitary filling, the ground state can be either a 
uniform Mott insulator with one particle per site or a 
charge ordered crystal with two particles on one sublat- 
tice and empty on the other one which is labeled as CDW 
H. There is a critical region around U ~ zV, where the 
two states with different translation symmetries become 
degenerate, and however, they are separated thermody- 
namically, i.e., any local perturbation cannot take one to 
the other. Correspondingly, the transition between them 
is a first-order one. It is noted that the aforementioned 



transition from the superfluid to CDW I state by tuning 
the chemical potential is of the weak first-order [21| . Far 
less attention has been paid to the region with zV ^ U 
by now, of which the details are given as part of the phase 
diagram in this article. To plot the ground state phase 
diagram, we focus on the case with small hopping and 
average particle density around or smaller than 1. For 
larger t or p, we expect no essentially new physics. 

This article is organized as follows. In Sec. [Hi we shall 
present the global phase diagram. The details of the 
order parameters will be discussed in Sec. IIIIl The con- 
clusion will be given in last section. 



II. GLOBAL PHASE DIAGRAM 

A. Classical case with t = 

We start from the classical case without hopping. The 
energy per site of ground state is a function of the particle 
numbers on the two sublattices, ua and n^, 

e'- >(nA,nB) = --[ua + riB) + —nAriB 

+ j[nA{nA - 1) + neins ~ 1)] , (2) 

where the coordination number z — 6 for the simple cubic 
lattice. The states can be labeled by {ua, ns)- The Mott 
states correspond to ua = ns , and the CDW states with 
nA ns break the translational symmetry, which is two- 
fold degenerate. In this article, we define the state (1,0) 
as CDW I state, and (2, 0) as CDW II state and we only 
consider ua > ns for the CDW states for convenience. 
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FIG. 1: (Color online) The phase diagram with zero hopping 
t = Q, where the states are labeled by the particle numbers on 
two sublattices (nA,ns) with the assumption ua >nB- The 
states with ua < ?^s can be obtained by inversion. The CDW 
states with one sublattice empty on the right side. Some 
states can only exist on the green thick solid line, e.g., the 
states (3,1) and (1,3). 

For /i < the ground state is a vacuum without any 
particles. As the chemical potential is increased (/i > 0), 
the particles are loaded into one sublattice to form a 
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charge ordered pattern with ua = I and = 0, i.e., 
CDW I state. If we further increase the chemical po- 
tential, more particles are loaded into the cubic lattice, 
which fill cither the empty sites if zV/U < 1 to form 
a uniform Mott state, or the occupied sites if zV/U > 1 
leading to a CDW II state. In the Mott state, each parti- 
cle interacts with its nearest neighbors, which effectively 
lowers the chemical potential to be /i — zV, and then the 
critical line between CDW I and Mott states is simply 
fi = zV. While that between CDW I and II states is 
a horizontal line ^ = U because the chemical potential 
only needs to compensate the on-site interaction U for 
adding new particles. 

Similarly, by studying the instability of adding parti- 
cle to a state {hatTib), one can determine all the phase 
boundaries between different classical insulating states, 
as shown in Fig. [T] There is a special vertical line 
zV/U — 1, on which many states can coexist with the 
same free energy. For example, on the boundary be- 
tween Mott and CDW II states, there are actually three 
macroscopic states, which are (1,1), (2,0) and (0,2). In 
fact, some of them only exist on this line in the ab- 
sence of hopping terms, e.g. states (3,1) and (1,3) on 
the boundary between Mott (2,2) and CDW (4,0) states. 
These degenerate states are macroscopic which cannot 
be transformed to each other smoothly by local pertur- 
bations, i.e., there are infinitely high barriers between 
these macroscopic states. 
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FIG. 2: (Color online) The phase diagram in the plane of 
^/U vs. zV/U for the extended Bose Hubbard model on a 
simple cubic lattice, where U = 40t and the lattice size is 
12 X 12 X 12. Four kinds of phases including the superfluid 
(SF), supersolid (SS), Mott and CDW states are identified. 
The CDW states are further classified as I and II by the filling 
numbers. The solid lines with circles are the phase boundaries 
calculated by QMC simulations, and the dotted lines are from 
the perturbation calculations. The green lines are the second- 
order phase boundaries, and the red ones are of the first-order. 
The inset is a zoom near the tip of CDW I lobe, where one 
can find a narrow region below the lobe corresponding to a 
supersolid phase with the filling less than one half. For details 
of order parameters one may refer to Figs.O [6] and [T] 



B. Case for finite t 

For the case with a finite hopping i, the particles 
(holes) adding to an insulating state can gain a kinetic 
energy, which results in the shrinking of insulating ar- 
eas in the phase diagram (Fig. ^ comparing with the 
classical case (Fig. [T|). In three dimensions, these mobile 
bosons condense at low temperature leading to the su- 
perfluidity, which enriches the phase diagram. There are 
four phases for the EBH model characterized by the fol- 
lowing three quantities: the particle density p, superfluid 
density ps, and static structure factor 5ff at momentum 

= (tt, tt, tt) with 
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where Wx,y,z are the winding numbers along x,y and z 
directions, /3 is the inverse temperature, and m = 2/t is 
the effective mass of the bosons. In the insulating Mott 
and CDW states, particles are localized by the interac- 
tion and the local particle number is quantized as in the 
classical case. The pure superfluid state has nonzero su- 
perfluid density ps and a vanishing static structure factor 
5^, while both are finite in the supersolid phase. 

Fig. [5] is the phase diagram determined by the QMC 
simulation with the worm update algorithm, where the 



solid lines with circles are the QMC results and the dot- 
ted lines are from the perturbation expansion in the 
strong coupling limit [l^ where the insulating states be- 
come unstable against adding or removing particles. It is 
seen that the perturbation results agree quite well with 
those of the QMC simulation in part of the phase dia- 
gram, but it is still not applicable in some regions since 
it cannot deal with the superfluid order. 

Comparing with the classical case in Fig. [1] the CDW 
I state is detached from its insulating neighbors, i.e., the 
Mott, CDW II and vacuum states. The upper boundary 
of vacuum state is actually lowered below p. = Q due to 
the hopping of bosons, which is not shown in Fig.[TJ The 
gaps between different insulating states are filled with the 
superfluid and supersolid phases. The lower boundary 
of CDW I state is a critical line on which there occurs 
a phase separation between the superfluid and CDW I 
phases, which breaks the U(l) gauge symmetry and the 
translational symmetry, respectively. The transition be- 
tween them is of the weak first-order across which the 
particle density and superfluid density have a jump. Con- 
sidering the correspondence between the EBH model and 
spin models, this transition is similar to the spin-flopping 
process in the two- or three-dimensional anisotropic XXZ 
model in the presence of magnetic fleld pointing in 2;-axis, 
which is equivalent to the EBH model in the hard-core 
limit [33l435| . As explained in Ref. [5|, this flrst-order 
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phase transition with a particle number jump is due to 
the fact that the CDW I phase is unstable toward the 
domain wall formation if the filling number exceeds L^^ 
though it is still stable against doping one hole. 

Doping particles upon the CDW I phase by increas- 
ing the chemical potential does not lead to a first-order 
transition as in the hard-core EBH model (or equiva- 
lently the XXZ model), where the particle-hole symme- 
try makes the upper and lower boundaries of CDW I 
phase identical. In case of the soft-core bosons, these ad- 
ditional particles can move upon the alternating charge 
ordered background with the effective hopping amplitude 
t^, which can Bose condense at zero temperature without 
destroying the staggered density order, and thus leads to 
a supersolid state with p > 0.5. This transition is of a 
second-order, as shown by a green solid line in Fig. [2] 
where all the second-order phase boundaries are in green 
that are distinct from the first-order ones which are col- 
ored in red. 

Continuously increasing the chemical potential upon 
supersolid phase, two different situations occur. (1) For 
zV < U, the particles like to occupy the empty sites 
which weakens the CDW order accompanied by the oc- 
currence of superfluid order. Until some critical filling 
p < I, the CDW order is completely destroyed and a pure 
superfluid state appears with the translational symmetry 
restored. The transition is of the second-order. (2) For 
zV > U, the additional particles are added to the occu- 
pied sites so that the CDW order is actually enhanced un- 
til finally entering into another insulating state, CDW II, 
through a first-oidei phase transition, of which the rea- 
son is the same as that from the superfluid with p < 0.5 
to CDW I. In this sense, the staggered density order in 
supersolid phase is inherited from the CDW I state, not 
related to the CDW II state, since doping holes to the 
CDW II state cannot smoothly result in the supersolid 
state. 



0.5 
0.4 
0.3 
0.2 
0.1 



Ps 



----A--- 



2|s.-^-- A----A A A-- 

>-^---G -e -0-- 



zV/U=0.3635 
H/U=0.145 - 



0.05 



0.07 



0.09 0.11 
1/L 



0.13 



FIG. 3: (Color online) The extrapolation of the order param- 
eters in the supersolid phase with p < 0.5. In the thermody- 
namic limit, the staggered structure factor S-w = 0.14, the su- 
perfluid density ps = 0.17 and the particle density p = 0.496, 
which indicates a supersolid state less than the half filling. 

As shown in the previous study [H, doping holes into 
the CDW I state may not lead to supersolid state with 



p < 0.5 on a two-dimensional square lattice, contrary to 
the case of doping particles. The reason is that the CDW 
I state has two-fold degeneracy. As long as enough par- 
ticles are removed, the insulating state becomes unstable 
towards the formation of domain walls between the two 
degenerate states, but can still be stable against losing 
one particle. This explains the particle density jump 
across the first-order phase boundary between the super- 
fluid and CDW I phases and that between the supersolid 
and CDW II phases. However this argument is invalid 
for V close to t that is around the tip of the CDW I lobe, 
where we show below that the kinetic energy gain by dop- 
ing one hole into the CDW I state is also linear in t, that 
can cause instability of CDW I state towards a supersolid 
state without the formation of domain walls. Suppose 
that one particle is taken from a CDW I state, the hole 
leaving behind moves in an effective staggered potential, 
roughly speaking, in one sublattice and (z — 2)V in the 
other. Solving this single particle problem, one estimates 
the kinetic energy gain AK is 



AK 



-tz 
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which is about — 2.88t in the cubic lattice near the tip 
of CDW I lobe where V ^ 2At (see Fig. [2]). As a con- 
sequence, a supersolid phase with p < 0.5 occurs, whose 
boundary is plotted in the inset of Fig. [2] To confirm, 
we also extrapolate the static structural factor, super- 
fluid order and particle density at point zV = 0.3635J7 
and fi — 0.145 to the thermodynamic limit (see Fig. [3]), 
which indicates a supersolid state with the particle den- 
sity less than half filling. For more information, one can 
refer to next section where we shall examine the order 
parameters in details. 

In the presence of a kinetic term, the vertical boundary 
at zV ~ U between Mott and CDW II states does not yet 
split but moves slightly to the CDW side. On this bound- 
ary, the free energies of both states are equal. However, 
they have different symmetries and are separated from 
each other thermodynamically. The transition between 
these two insulating states is of the first-order, which 
is similar to the conventional liquid-solid phase transi- 
tion. Further doping particles into the Mott and CDW 
II states leads to the reentrance of superfluid and super- 
solid states, respectively, with density p > 1. 

Because the particles or holes can hop between near- 
est neighbors on a Mott background, which leads to a 
kinetic energy gain linear in t, however they can only 
hop among next nearest neighbors in a staggered CDW 
background that gives rise to a hopping energy propor- 
tional to t'^/V. Then, as t increases the Mott region 
shrinks much faster than the CDW II region, which re- 
sults in the mismatch of phase boundary between them 
as zV ^ U{see Fig. [2]). Since doping the Mott state with 
holes leads to a superfluid state, one expects a boundary 
between the superfluid state and the CDW II state in the 
critical region at the lower end of the vertical boundary. 
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This phase boundary is again of the first-order. At the 
upper end, the extension of the vertical boundary sepa- 
rates the superfluid and supersohd phases in the way of 
the second-order phase transition. 



III. ORDER PARAMETERS 
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FIG. 4: (Color online) The average particle density p as a 
function of chemical potential ^ for different V and U — AOt. 
The first plateau at p — 0.5 corresponds to the CDW I phase, 
and the second one corresponds to the Mott phase for V/t = 
3, 4, 5 and the CDW II phase for V/t — 9. The system size is 
12 X 12 X 12. 

In this section we give the details of the parameters p, 
Ps and S-g for different p and V. Fig. S] is the density 
profile as we vary the chemical potential for several fixed 
values of V , where the plateaus correspond to the incom- 
pressible states, i.e., the Mott, CDW I and II states, with 
vanishing isothermal compressibility kt = p^^dp/dp. 
Note that these insulating states correspond to two sin- 
gle points p = 0.5 and p = 1 in Fig. [5l which is the plot 
of superfluid density ps and as functions of particle 
density p. 

For p < 0.5, the particles Bose condense to form a 
pure superfluid state with a nonzero ps but vanishing 
as shown in Fig. [SJ When the particle density reaches a 
commensurate value p = 0.5, a plateau appears implying 
Kt — which corresponds to the incompressible CDW 
I state with translational invariance broken. This tran- 
sition is of the first-order since the particle density, as 
a first-order derivative of free energy with respect to the 
chemical potential, is discontinuous. This is also reflected 
in Fig. [5] as the fact that a segment of particle density 
below the half filling is inaccessible, and jumps at 
p = 0.5 to a finite value from zero. As the particle den- 
sity exceeds 0.5, kt becomes a finite positive value again. 
At the same time the superfluidity appears continuously 
upon the CDW I state to form a supersolid. The cor- 
responding transition is of the second-order, since the 
particle density p is continuous but kt jumps from zero 
to a finite value. 

Between the two plateaus at p = 0.5 and p = 1, the 
slope of p{p) curve for V/t = 3,4 and 5, i.e. zV < U, 
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FIG. 5: (Color online) Superfluid density ps and staggered 
structure factor Sjf as functions of the particle density p at 
U/t = 40. Note that the plateaus of constant p in Fig. |4]are 
only single points here. 

shows another jump, implying that kt is discontinuous 
and a second-order phase transition from the supersolid 
to superfluid phase occurs. This transition is manifested 
in Fig. [5l where ps keeps finite in the whose region 
0.5 < p < 1.0, but Sa vanishes at some critical values in 
between. Further increasing the chemical potential, the 
system enters into the Mott phase, which corresponds to 
the plateau with p = 1 in Fig. 21 through a second-order 
phase transition. At the transition point ps vanishes as 
shown in Fig. [5] 

For V — 9t, i.e., zV/U > 1, the second plateaus in 
Fig. m corresponds to the CDW II state, into which the 
system enters directly from the supersolid phase as p is 
increased. The transition is a first-order one as shown in 
Fig. 13] as the jump of particle density. It is also reflected 
in Fig. [5] where a segment of particle density p is not 
accessible before p reaches 1. 

In Fig. ini we plot the normalized order parameters 
Ps = Ps/Ps,max and = Sjf/Sjf^max as functious of 
zV/U around the CDW I lobe, where Ps,max ~ 0.6 and 
Sjf^max ~ 0.2. In the plot, we take three characteristic 
values of the chemical potential p/U = 0.160, 0.155 and 
0.145 from the top to the bottom panel. In the shaded ar- 
eas, when V increases, the superfluid density ps decreases 
while the static structure factor 5^ increases, and both 
Ps and Sj; are nonzero in this region indicating a su- 
persolid state. The supersolid area becomes more and 
more narrow as p decreases, until finally shrinks to a 
point in the zV/U axis when p/U 0.140, which im- 
plies that two second-order phase boundaries merge into 
one first-order phase boundary across which the order 
parameters ps and 5*^ are both discontinuous. The cor- 
responding density profiles are plotted in Fig. [3 where 
we observe that the particle density is always smaller 
than 0.5 for p/U — 0.145 corresponding to the bottom 
panel in Fig. [S] until it enters into the CDW I phase, that 
implies the supersolid state can exist below the half fill- 
ing. For p/U = 0.160, the situation is different where 
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FIG. 7: The particle density vs. zV jJJ for /i/C/ = 0.160, 0.155 
and 0.145. 



0.359 0.36 0.361 0.362 0.363 0.364 0.365 
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FIG. 6: The normalized order parameters ps and as func- 
tions of zV jlJ . The upper, middle, and bottom panels cor- 
respond to [ijU = 0.160, 0.155 and 0.145, respectively. The 
supersolid phases occur in shaded areas. The lattice size is 
12 X 12 X 12. 



p decreases as V is increasing, which can be intuitively 
attributed to the loss of effective chemical potential due 
to nearest neighbor interaction V in the mean-field level. 
The case of pjU — 0.155 shows an intermediate behav- 
ior, for which the particle number first decreases, then 
increases slightly larger than 0.5 and finally reaches the 
CDW I state. 
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IV. CONCLUSION 



In this article, we present the global phase diagram 
(Fig. of the 3D extended Bose Hubbard model. The 
EBH model exhibits four kinds of ground states, includ- 
ing (1) the Mott state without breaking any symmetry, 
(2) the CDW I and II states with translational symme- 
try broken, (3) the superfiuid with U(l) gauge symme- 
try broken, and (4) the supersolid with both symmetries 
broken. By using the QMC simulation as well as other 
analytical tools, we identify the transition type between 
these phases. Among them, the first-order phase bound- 
ary includes those between the superfiuid and CDW (I 
and II) states, and between the Mott and CDW II states. 
The other boundaries are all continuous. The critical re- 
gions for zV jU ^ 1 and the tip of the CDW I lobe are 
examined in detail. We demonstrate that in the present 
3D EBH model, the supersolid phase with p < 0.5 can 
appear in a small region near the CDW I lobe where the 
hopping amplitude t is comparable to the nearest neigh- 
bor interaction V. In this region, the general "domain 
wall" argument for the nonexistence of supersolid state 
with p < 0.5 is no longer applicable, since it is based on 
the assumption oft-^V. 
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